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Abstract 

Using molecular dynamics simulations, we determine the linear and nonlinear viscoelastic 
properties of a model polymer melt in the unentangled regime. Several approaches are compared 
for the computation of linear moduli and viscosity, including Green-Kubo and non equilibrium 
molecular dynamics (NEMD). An alternative approach, based on the use of the Rouse modes, 
is also discussed. This approach could be used to assess local viscoelastic properties in inhomo- 
geneous systems. We also focus on the contributions of different interactions to the viscoelastic 
moduli and explain the microscopic mechanisms involved in the mechanical response of the melt 
to external sollicitation. 
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I. INTRODUCTION 



The response of polymer melts to mechanical perturbations, involving either oscillatory 
or steady flow, is of great practical importance, and has been the object of extensive 
experimental and theoretical studies. 3,13,0] This response is well known to be viscoelastic, 
i.e. the storage and loss moduli exhibit a strong frequency dependence, and nonlinear, 
with a typical shear thinning behavior for the viscosity. 

On the simulation side, the steady state viscosity and shear thinning effects have 
been studied quite extensively in various configurations for model systems. An extensive 
review of recent work is given in reference. A spectacular success was the obtention 
of the Rouse- Reptation (or unentangled-entangled^ crossover in the rheological behavior 
for model polymer melts of the bead spring type. 6] This crossover is obtained for chain 
lengths of the order of = 100 monomers. According to the popular wisdom in the field, 
melts with A^ < 100 should therefore be amenable to a description in terms of the Rouse 
model, which is considered as a reasonable phenomenological description of unentangled 
melts. 

Investigation of frequency dependent response are much less numerous than for steady 
state viscosity. In fact we are aware of only one recent study, with objectives quite sim- 
ilar to those of the present article. Our aim is to investigate the mechanical response of 
model polymer melts submitted to steady or oscillatory shear, in order to obtain a charac- 
terization in terms of frequency and amplitude of the solicitation. In view of the numerical 
cost of such calculations, our study will be limited to short chains, but will explore several 
values of amplitudes and frequencies, concentrating on relatively low frequencies. 

As it turns out, a direct assessment of mechanical properties using non equilibrium 
molecular dynamics (NEMD) is very costly from a computational viewpoint. Hence it 
is desirable to explore methods that could provide the same information with a lesser 
computational effort. We will in particular explore the possibility of obtaining viscoelastic 
properties directly from a study of Rouse modes. Those, being single chain properties, 
offer a much better statistical accuracy than the stress itself, which is a global property 
of the system. 
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The systems under study are briefly described in the next section. We then discuss 
the steady state viscosity, both in the hnear and nonhnear regime. We use three different 
methods to determine the viscosity - NEMD simulations, equilibrium Green-Kubo ap- 
proach and show how the viscosity can be obtained in a third way from the analysis of 
equilibrium Rouse modes, provided the contribution from short times is correctly taken 
into account. We then turn to the study of oscillatory strains by means of NEMD simula- 
tions. The conditions for linear response at low frequency are discussed, and the various 
contributions to stress response - storage and loss - are estimated. Again, the results are 
compared to an equilibrium analysis based on the Rouse model. 



II. SYSTEM DESCRIPTION AND METHODS 



The chains are modelled by an abstract and generic, though well studied, bead spring 
model - the rather common " Lennard- Jones + FENE" model. |1] All monomers in the 
system are interacting through the Lennard- Jones potential: 



r4£((a/r)i2_(a/rf), r<r, 
Uij{r) = < (1) 
10, r > Tr 



where = 2.5cr. Neighbor monomers in the same chain are linked by the FENE (Finite 
extension non-linear elastic) potential: 

UpENE{r) = ^RoH'^-{^)'), r<Ro (2) 

where Rq = 1.5a and k = SO.Oe/a^. The typical size of the studied systems was of 1280 
and 2560 particles (for chain lengths of 10 and 20) in a periodic cubic simulation box. 
The temperature of the melt was fixed in all simulations at ksT = 1 and the density 
was p = 0.85. For the non equilibrium runs we use the Lees Edwards periodic boundary 
conditions coupled with the SLLOD ^1 equations of motion. This algorithm allows us 
to simulate an unconfined periodic system under shear. The stress tensor in the melt 
is calculated as (Ja/sit) = —y '^r'^jF^-, where F^- the (3 component of the force between 
particles i and j. The kinetic contribution of the stress (oc '^vfv^) was also evaluated 
and found to have a negligible contribution to the stress as expected at high density. The 



simulation time step was taken to be 5t = 0.005tlj, where tlj is the Lennard- Jones time. 
All times in the resu 
the LAMMPS code. 



ts are in units of LJ time. The simulations were performed using 
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III. STEADY STATE SHEAR VISCOSITY 



NEMD determination 



The most direct method to obtain shear viscosities is undoubtedly the nonequilibrium 
approach that generates homogeneous, planar Couette flow using Lees-Edwards boundary- 
conditions and the Sllod algorithm. For a shear rate 7 and a velocity profile Vx{y) = jy, 
the viscosity is given by 

(3) 

7 

Where a{t) is the xy component of the stress tensor in the sample. The viscosity obtained 
by this method depends on the shear rate, generally decreasing with 7 (shear thinning). 
An extrapolation is required to estimate the zero shear rate value. Shear thinning gener- 
ally takes place when 7 > 1/tc, where Tc is a characteristic relaxation time of the polymer 
melt (usually, for unentangled melts, the Rouse time t^, as can be seen in fig. [T]). 

Precise measurements for low shear rates are very time consuming as substantial statis- 
tics are needed for the accurate determination of (cr(t)) which has a small value for small 
shear rates. As the relaxation time of the polymer chains increases rapidly with the chain 
length (tc oc N^), for long chains very low shear rates 7 < I/t^, should be used to reach 
the Newtonian plateau in ^7(7). This is illustrated by the large error bars on the low shear 
rate side of figure For chains of lengths = 10 and = 20, our extrapolation at zero 
shear rate is consistent with the scaling expected for unentangled melts, rj (x N, and with 
earlier estimates found in the literature. 
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FIG. 1: Shear viscosity as a function of shear rate 7 measured in simulations of planar flow 
(equation |3J . Shear thinning takes place for approximately 7 > 0.01 for = 10 and 7 > 0.0025 
for N = 20. The estimated Rouse relaxation times for the two chain lengths are respectively 
about 100 and 400. Stress was measured every 2tlj for a time ranging from ISOOOr^^j for high 
shear rates to lOOOOOr^j for low shear rates. Error bars are estimated using block averages for 
data correlation. jin| 

B. Green-Kubo approach 

The problem of extrapolating the shear rate dependent viscosity value does not occur 
in equilibrium methods. The zero shear rate viscosity is given by the integral of the stress 
correlation function (Green Kubo relation) 

V 



Git) 
V 



G{t)dt 



(4) 
(5) 



with T and V the melt temperature and volume. Obtaining an accurate value for the 
slowly decaying stress correlation function involves very long runs in order to have suffi- 
cient statistics. Assuming for simplicity that a{t) obeys Gaussian statistics, the standard 
error in the time correlation function at long times is given by 



5((a(t)a(0))) 



(6) 
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where Tc is a characteristic correlation time of the data. Given that the stress correlation 
is a rapidly decaying function (see fig. HI (a^) ^ 10^(cr(10)cr(0))), this means that if 
we want to estimate correlations for a correlation time ^ 100 it would be necessary 
to run a simulation of 10^ (units of LJ time) in order to obtain a relative precision of 
~ 10% for (o"(100)o"(0)). The situation in the polymeric case is particularly unfavorable, 
since the large viscosity is obtained as the product of a relatively low modulus and a long 
relaxation time. Hence the stress correlation function is long ranged in time, but has a 
low amplitude, easily masked by noise associated with the rapid pair-potential part of 
the stress that does not involve polymer chain relaxation. This estimate is indeed very 
pessimistic and should be seen as an upper bound of the possible error in the Green Kubo 



formula. [ll| In practice the stress autocorrelation function has several contributions with 
different weights and correlation times so that the exact error in the Green Kubo integral 
is difficult to evaluate and depends on the chain length. In our systems of relatively short 
chains the uncertainty of the Green-Kubo formula concerning viscosity is of the order of 
30% when the integral is carried out to several relaxation times and grows larger as the 
correlation function is further integrated (see figure Ej). 

We did not examine the infiuence of system size on the statistical accuracy of our 
data. If a larger system should diminish stress fiuctuations (by a factor oc V^, where 
is the number of particles), a similar effect is expected by longer simulation time (a factor 
oc x/T). As the computational effort for larger system increases roughly as A^logA^, we 
do not expect to gain better precision for less computational time this way, so we did not 
examine the trade off between system size and run time. 

We are aware of several determinations of polymer melt viscosity using the Green-Kubo 



approach. In 



12j | a good agreement with NEMD results was obtained for short simulation 



times without discussion on uncertainty. Another calculation was made in [15], where the 
authors reported the existence of a large noise in the correlation function. This problem 
was solved by performing running averages to smooth the data, a method that can reduce 
the noise due to rapid bond vibrations but whose effect on the the intrinsic statistical 
accuracy of the stress correlation is not obvious. Our data still indicate a large error bar 
(30%), when this error bar is estimated from the three independent components of the 
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stress tensor. A viscosity determination based on a Green-Kubo formula in terms of an 
Einstein relation was presented in unfortunately with relatively little details that 
would allow us to compare with our results in terms of efficiency and accuracy. It is clear 
that the Green-Kubo formula remains the only exact way of determining the viscosity 
from equilibrium simulations, and should be used whenever an "exact" result is required. 
Although this is feasible with a large computational effort. a detailed report on its 
accuracy for polymers is still missing (see however ^11]). Therefore it seems interesting 
to discuss an alternative, faster, approach that can be used for example in comparative 
studies at a moderate computational cost. 



C. Rouse modes 



It is well known that, at least at a qualitative level, the Rouse model can account for 
the viscoelastic behavior of unentangled polymer melts. Hence, it is tempting to attempt 
to bj^ass the difficulty in obtaining the viscoelastic properties stricto-sensu by directly 
using this model. The large uncertainties discussed above are highly reduced, when we 
turn to the calculation of single-particle (or single-chain) correlation functions, which are 
the essential ingredient of the Rouse model. As the final result is an average over M 
separate functions the standard error at long times should be ~ {2tc/ MtrunY^'^ ■ 

The spirit of the Rouse model consists in assuming that the melt mechanical behavior 
is dictated by the relaxation of a single polymer chain, the influence of inter chain in- 
teractions being limited to the phenomeno logical friction constant. Consistent with this 
assumption, the mechanical stress can be calculated from the Rouse modes of the chains 

where 

^.W = ^E-"W-^(^^^^^#^)' P = 0,...,Ar-l (8) 

n=l ^ ^ 

are the Rouse modes, with the chain length, p the monomer number density and r„(t) 
- the position of the n-th monomer in the chain at the time t. Assuming, again in the 
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spirit of the Rouse model, independent Rouse modes, the stress correlation from equation 
can be rewritten as a function of the equilibrium correlation functions for individual 
Rouse at equilibrium: 

G{t) = -i^(a.,(t)a.,(0)) (9) 



y 1 r°"" fpkBTy^{Xp,it + r)X,y{t + T)) (X,,.(r)X,,(r)) 



^ l_ (pknTy 1 ^^ X^^,it + T)X;yit + T)X^',ir)X^'M 

_ J_ r"", pksT 1 ^""^ x;,{t + T)x;y{t + T)x;,{T)x;y{T) 

_ p/c^T (X,,.(t)X,,(t)X,,.(0)X,,(0)) 

The viscosity is then calculated using equation With this method we observe a 

substantial gain in precision (see figure |2I) and the values obtained are smaller than the 
non equilibrium estimates, by about the same amount for the two different chain lengths 

{VNEMD -Veq^ 5). 

This result is not surprising given the simplifications of the latter calculation. The 
interactions between chains in the melt are not taken into account and the Rouse model 
cannot yield information about stress relaxation on very short "non - polymer" time 
scales. In order to get a better understanding of the deficiencies in the calculation using 
the Rouse model, we first compute the individual relaxation times of the modes. The 
relaxation times, shown in figure El are extracted from the exponential decay of these 
correlation functions. The Rouse scaling Tp oc l/p^ is well obeyed for the first modes. 13 1 
While this is a typical result for the bead-spring model, we note that using more 
detailed atomistic simulations larger deviations from the Rouse scaling can be observed, 
especially in higher modes, jlj] For both chain lengths, the relaxation time of the fastest 
mode was found to be Tjy^i ~ 2. 

It is clear that the Rouse calculation cannot account for the contribution to the viscosity 
associated with time scales shorter than r^r.!. On such short time scales the Green-Kubo 
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FIG. 2: Integrals of the stress correlation function G{t) for = 10 obtained from the global 
stress (Green-Kubo) and the Rouse modes stress. A distinct plateau that stays stable for very 
long times (~ 50t/j) is clearly reached with the Rouse modes formulation. For comparison the 
Green-Kubo integral value has important fluctuations for times longer than several times the 
Rouse relaxation time as predicted by the uncertainty calculation. Global stress was measured 
every O.Olr^j for a time span of 25000rij and G{t) was calculated by taking every point as a 
starting state and averaging. Such fine sampling is needed to capture the fast oscillations in 
G{t) shown in fig. |3J Rouse modes were measured every CSr^j for every chain for a time of 
15000rij. Than G^'"""'{t) was calculated averaging over starting states and chains. 

viscosity integrand can however be obtained with high accuracy, as illustrated inEJ This 
stress-stress correlation function first decreases rapidly, and displays a short time damped 
oscillatory behavior, with a time constant smaller than 0.5. Similar stress oscillations were 
reported for n-alkanes as well as for a bead-spring melt, jisl. Q] In our case, a detailed 
study of the FENE bonds and Lennard- Jones forces contributions to the stress, shows 
that these oscillations are due bond vibrations, and that their frequency is close to the 
intrinsic frequency of the FENE bonds. 

Integrating this short time stress correlation function gives the contribution to the melt 
viscosity of the rapidly decaying part of the global stress, unaccounted for in the Rouse 
model. 




(14) 



(15) 
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FIG. 3: Relaxation times of the Rouse modes of the chains versus the mode number. Relaxation 
times were estimated by an exponential fit of the normalized correlation function of the mode, 
which had in all cases an exponential behavior. Dashed lines show the Rouse theory prediction 
Tp oc followed closely by the first several modes as discussed in 

Adding the contribution of the second term of equation to the viscosity from the 
Rouse model provides a significantly more accurate estimate of the viscosity, compared 
to that calculated from non equilibrium runs (table The short time behavior of the 
stress correlation function is the same for chain lengths = 10 and = 20, which is an 
evidence that this short times contribution is independent of chain length and represents 
the missing "non-polymeric" part of the mechanical properties of the melt. It provides a 
contribution to viscosity of the order of the viscosity of a simple L J fluid at the considered 
density and temperature. 

Being chain length independent, short times stress relaxation represents a part in 
viscosity that diminishes with increasing chain length, about 40% for A^ = 10 and ~ 30% 
for A^ = 20. Still it cannot be neglected for predicting mechanical properties of a melt of 
unentangled chains. It is equally important for long chains at large shear rates, when the 
melt viscosity becomes comparable to the one of a simple fluid and can explain deviations 
from the stress-optic rule.j?! 

The presented method provides a way to rapidly obtain an estimate of the viscosity as 
it requires much less computational time than an accurate Green Kubo measurement. It is 
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FIG. 4: Global stress correlation functions for short times. Correlations of the pair Lennard- 
Jones and FENE bond interactions contribution to the global stress are also shown. Dashed 
line shows a fit of the global stress correlation to indicate the relevant decay time scales and 
oscillation frequency. 



Viscosity 


iV = 10 


N = 20 


Green-Kubo estimate 


8±4 


N/A 


Rouse modes value 


4.88 ±0.28 


10.2 ± 1.0 


Short times correction 


3.38 


3.43 


NEMD extrapolation 


8.2 ± 1.3 


14.9 ± 1.9 


Corrected Rouse value 


8.26 ±0.28 


13.6 ± 1.0 



TABLE I: Viscosity values from the different methods. Uncertainties in Green-Kubo and 
Rouse modes G{t) integrals are estimated by evaluating the {axy{t)axy{0)) , {crxz{t)crxz{0)) and 
{cryzit)cryz{0)) integrals. 

still less straightforward and one needs to monitor (where is the chain length) Rouse 
modes for every chain to build the correlation function. The method's main advantage 
is that, being based on a single chain quantity, it is local in nature and can be used to 
assess local viscoelastic properties in inhomogeneous systems such as confined melts or 
melts with filler particles. 



Full stress correlation : 
47exp(-t/0.1)cos(43t) + 21exp(-t/0.1) : 
Bond stress correlation : 
Lennard- Jones stress correlation : 
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IV. ELASTIC MODULI 



A. Method 

Like the viscosity, elastic moduli can be obtained using either equilibrium or nonequi- 
librium simulation. If one is interested only in linear response properties, the moduli can 
be obtained by Fourier transforming the Green-Kubo integrand G{t) (equation in the 
form 

POO 

G'{uj) = LU dtG{t)smu;t (16) 
Jo 

POO 

G"{uj) = uj I dtG{t)cosujt (17) 







This equilibrium determination, however, suffers from the same drawbacks as the Green- 
Kubo determination of the viscosity, i.e. a long time, small amplitude tail of G{t) has 
to be known very accurately to obtain reasonable results. That is why we did not show 
any result concerning elastic moduli issued by a Green-Kubo method. Instead, the only 
practical way of using equations II (il and fTTI is to start from a "model" calculation of G{t), 
in the sense of the Rouse modeling described in the previous section. 

To obtain information beyond the linear regime, the only possibility is to effectively 
do NEMD and submit the sample to an oscillatory strain, using the standard SLLOD j9| 
algorithm. The strain is given by 

7o(^) = 7oi^sinciJt =^ 7o(^) = —70 cos cut (18) 

And we define the response of the system by the usual formulae for the stress a{t) 

cx{t) = [ G{t-t';-fo,uj)io{t')dt' (19) 



A dependence on the shear amplitude and frequency is indicated in the response function 
G, to recall the possible existence of nonlinear effects. The frequency dependent moduli 
are defined from the Fourier component of a{t) at the imposed frequency u 

cr(t) = 7o(G"(ti;, 7o) sincut -|- G'"(co', 7o) costut) + harmonics at 2uj,3uj.. (20) 
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The moduli are formally given by the Fourier transforms of the response function 

POO 

G'(a;;7o) — u dtG{t; ^o,<^) sin ut (21) 
Jo 

/•oo 

G"{cu;^o) = cu dtG{t; ^0,1^) cos cut (22) 
Jo 

In practice, G' (resp. G") is extracted from the time series for the stress by multiplying 
the signal by cos(a;t) (resp. sin(a;t)), i.e. 

G'(ci;,7o) / dtcos^ujt = / dta{t) cos ut + G" {uj) I dtcosutsinut (23) 

Jo 7o Jo Jo 

with Tr the length of the simulation run. Thus we obtain the storage and loss moduli as 
a function of stress : 

G'(a;,7o) = (-^ ^ dta{t) cos ut + G''(u;;^of^^) (24) 

G'\u) = ^— 1^ (l£^rfta(t)sina;t + G'(a;;7o)^^^) (25) 

In the above formulae potential harmonic terms were ignored for simplicity, their contri- 
bution being of order l/T^. In the limit ^ l/o; one has simply 

2 /"^ 

G"(a;;7o) = / dta{t) cos ut (26) 

TrJO Jo 

G"{u;-fo) ^ 7^ I dta{t)sinut (27) 

In order to elucidate the role of the different interactions for the elastic moduh we focus 
on their respective contribution. The stress can very generally be separated into an 
intramolecular stress component associated with FENE bonds and intra chain Lennard- 
Jones forces, and a intermolecular stress component associated with inter-chain Lennard- 
Jones interactions. In the following the moduli issued from NEMD simulations will be 
discussed in terms of these two separate intra and inter molecular contributions. During a 
NEMD run the full stress in the system as well as its inter and intra molecular components 
are stored and then elastic moduli are calculated as described above. 



13 



B. Results 



1. NEMD Results 

The first question that we investigated is the extent of the hnear regime, in terms of 
the strain amplitude. Nonlinear effects can in principle be detected by a dependence of 
G{uj) on 7o, or by the presence of higher harmonics in the stress signal. 

There is a clear softening of the response at frequency uo as amplitude is increased. 
This softening is obtained above a value of the strain rate 7oCt; of the order of at low 
frequencies, as illustrated in figure The situation is very similar to the shear thinning 
behavior of the viscosity, namely the relevant parameter is the shear rate rather than the 
strain amplitude or frequency. 
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FIG. 5: Dependence of elastic moduli on strain rate amplitude 70a; for = 10. Elastic moduli 
are normalized by their lowest shear rate amplitude value. For all frequencies the linear regime 
extends at least to 70W ~ 1/t/j, for uj ^ I/t/j the linear behavior breaks down for ^qlo oc oj. The 
same behavior was measured for = 20 (data not shown) 

For frequencies uo ^ I/t/j, the softening is observed for values of ''Jquo significantly larger 
than I/t/j. The behavior of chains in this frequency range is further discussed below. 
Harmonics were detected in the time series for a{t) at high values of the strain amplitude 
7o > 2.5 at a frequency of Stu, where uj is the solicitation frequency, for all frequencies. For 
symmetry reasons the stress must be an odd function of the strain, so that the response 
function G is an even function of 7. Hence harmonic contributions are observed only 
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for odd multiples of the solicitation frequency. The amplitudes of the harmonic terms 
in the power spectrum are small, about 8% and 3% of the uo peak for respectively G'{uj) 
and G"{uj). In the molecular stress harmonics are more visible with about 17% and 
5% for the storage and loss moduli, respectively. We have not been able to distinguish 
harmonics for strain amplitudes 70 below 2.5. The observation of harmonics can thus 
be attributed to physical extension of the chains in which the non linear terms in the 
interaction potentials become inevitably important. Given this preliminary investigation 
of the non linear regime, we choose first to explore the linear response and place ourselves 
at shear rates 700^ < 1/tr. 

Figures El and [7| display the frequency dependence of the elastic moduli. As discussed 
above, the important uncertainty in G{t) obtained from equilibrium calculations imposes 
the use of non equilibrium methods. 

0,1 



0,05 

12 3 4 5 6 

FIG. 6: G'{uj) for N = 10, measured by NEMD. Contributions from the intra-molecular and 
inter-molecular forces are shown. Intra molecular forces are important for stress storage up to 
high frequencies. As > l/xp, the modes Xj, i < p are "frozen" and behave like stiff springs 
and store stress efficiently, so that the intra-molecular component of G'{io) grows with frequency. 

The contributions of the different interactions to the elastic moduli (calculated by 
eqn. (j26p and (|27p) can be examined by measuring the different contributions to the 
global stress. The loss modulus, which is a measure of stress dissipation in the melt has 
several contributions depending on the time scale (fig. [7|): at short times stress is relaxed 
through the pair interactions between monomers, in a "liquid like" manner. At longer 
times, relaxation of chain fragments of length Np = 1,2...N, take place on increasingly 
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FIG. 7: G"{lo) for N = 10, measured by NEMD. Contributions to the moduli from the intra- 
molecular and inter-molecular stress components are also shown. For to > I/tr the inter- 
molecular contribution to the loss modulus grows larger than the molecular component and the 
system crosses over to a liquid-like regime. 

larger time scales. Knowing that a mode p can be viewed as the relaxation of a subchain 
of Np = N/p monomers, for a given mode p relaxing over Tp, if the frequency is such 
that uj > 1/Tp, the mode cannot relax over one oscillation and does not take part in the 
stress relaxation, meaning that stress is relaxed on scales smaller than N/p monomers. 
This leads to the Rouse model prediction that, if chain relaxation was the only process 
in the melt, the loss modulus had to decrease at high frequencies. This decrease was 
not observed in our model melt (fig. Hj), due to non-polymer relaxation. As non-polymer 
relaxation we refer to the relaxation of stress occuring on a time scale smaller than the 
relaxation of the fastest Rouse mode and independent of chain length. Results show 
that at high frequencies the behavior of the loss modulus is dictated by inter molecular 
interactions, the inter molecular stress component is dominant. There is a crossover from 
polymer melt behavior where stress dissipation is carried out mainly by chain relaxation 
to a behavior of a liquid of interacting chains that do not have time to significantly change 
their conformation over one period. In this regime there is no stress dissipation due to 
internal polymer chain relaxation and the increase of the loss modulus with frequency is 
entirely due to Lennard Jones pair interactions on very short time scales. The situation is 
somewhat different for the storage modulus, where simulations show that internal polymer 
chain interactions and chain modes are important for the elastic response of the melt up 
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to high frequencies (fig. For a simple hquid the value of G'{u!) is very small in the 
considered frequency range, so it is not surprising that its value for the melt is due to the 
chains. The stress storage thus takes place to a large extent in the slowly varying chain 
conformations. Over a large frequency range the slow vibration modes act as an energy 
reservoir and the higher the frequency, the more the chains remain rigid at the time scale 
of a single period and thus cause the increase in the storage modulus with frequency. The 
melt exhibits an elastic behavior due to the chains that grows stronger with frequency, 
in fact the mechanism exposed for the loss modulus can be applied the other way around 
for G'{uj). A given mode p goes rigid as u > 1/rp and thus stores stress (rigid behavior) 
instead of relaxing it (liquid behavior). As > l/r^, the modes Xj, i < p are "frozen" 
and take an important part in energy storage. For very high frequencies u > 1/r/v-i the 
chain behaves more like a spring than a fiexible polymer. 
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CO 

FIG. 8: Storage modulus from NEMD and equilibrium Rouse modes calculation for chains of 
length = 10 and N = 20. For high frequencies the response of both systems: = 10 and 
= 20, becomes identical. 
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2. Calculation from Rouse modes and comparison to NEMD data 

An analytical calculation of G'(uj) and G"(uj) can be done using the Rouse model 

nn 

[2, y| . For the bead-spring polymer melt studied here these functions can be estimated 
using G^"^^'^{t) (eqn. via the integrals in equations fIT\i and ()22|1 . Given that, as 

17 



1 



0,01 



0,1 




0,001 



0,001 0,01 
CO 



0,1 0,001 0,01 



0,1 



FIG. 9: Elastic loss modulus from NEMD full stress (left, = 20 - squares, = 10 - circles), 
NEMD intra molecular stress (right, = 20 - squares, = 10 - circles), corrected equilibrium 
Rouse modes calculation (eqn. I32|) (left, A^ = 20 - stars, A^ = 10 - small circles) and Rouse 
calculation without short times corrections (right, A^ = 20 - stars, A^ = 10 - small circles). The 
uncorrected Rouse calculation fits well the intra molecular moduli component. Stress relax- 
ation for high frequencies is identical for the two systems, dictated by short subchain and inter 
molecular forces relaxation on short time scales as discussed. 

discussed above, the stress storage is dictated by slow, "frozen" chain vibration modes for 
the whole frequency range, the melt elastic response is well reproduced by this calculation 
(see fig. IHl). As shown in fig. IHl the Rouse approach leads to an underestimate of the 
loss modulus, especially at high frequencies {uj > 1/t/j). The increase in G"{uj) for high 
frequencies is due to short time non-polymeric stress relaxation discussed in the previous 
section. As discussed in the first part concerning viscosity, this contribution cannot be 
predicted by G^°^'^^{t) that takes into account only chain vibration modes. As no inter 
molecular forces whatsoever can be taken into account by the single chain Rouse model, 
QRouse^^^ provides a good approximation of the loss modulus calculated from the intra 
molecular stress component and shows a discrepancy with the full stress loss modulus 
that grows with frequency. In order to obtain an equilibrium single chain estimate of the 
mechanical behavior of the melt for all frequencies, we estimate the short time corrections 
to G^"^^'^{t) in a manner similar to the one used for viscosity. We use the fit of the short 
times stress correlation function (fig. 0]) to calculate the short times correction to the 



18 



elastic moduli in the frequency domain by equations 1)2111 and ()22j) . Writing the short 
time stress correlation in the form 



Gf^st(^t) = Ae-^l^' cos + Ee"*/"^ (28) 



leads to 



oo 



G'{uj) = uj / G^°'''\t)suiujtdt (29) 







+ 77 , , /. . , o^'^-2 + T-^. FTn^ + ^ -1 , .9-2 (^0) 



2 V1 + (^ + ^)'t2 1 + (o;-r])2^2y l + t^2^2 

/>CXD 

G"(cj) = UJ G^°'''\t) cos ujtdt (31) 







^ ^ Vl + (^ + ^)V ^ l + (a;-fi)2riV +^l + a;2r| ^^^^ 
where we determine the parameters A, B, Q, ti and T2 from the stress correlation function 
{A = 45, n = 43, n = T2 = 0.1 and = 21 for both = 10 and = 20). Adding 
these terms to the equilibrium Rouse modes loss modulus gives a much better estimate 
of G"{u!), producing the curves referred to as corrected Rouse (fig. El corrected Rouse). 
The correction concerning G'{u!) is negligible for the frequency range studied here and 
the corrected curve falls on top of the original Rouse curve. This is the expected result 
knowing that, as already mentioned, the liquid like interactions that dominate at short 
times participate in stress storage only at very high frequencies. We find the expected 
linear dependence of G"{uj) for < l/r/j (fig. ITT|l . the slope being within error bars 
the value of the viscosity estimated by planar Couette flow simulations. The storage 
modulus has, as expected, ~ u"^ behavior at low frequencies (cu < 1/t^) (fig. ITUI) . The 
Rouse theory predicts a cross over towards a oc ^/uJ behavior for higher frequencies 
Our simulations show that G'{u!), estimated by NEMD and Rouse modes measurements, 
grows slightly faster than ^/uJ at high frequencies (fig. fTIH) . We can relate this to the 
discrepancy between theoretical and measured modes relaxation times and argue that the 
"mean field" presence of multiple chains in our vibration modes determination promotes 
more efficient stress storage in the melt at high frequencies. Simulations show that G"{u!) 
does not follow the a/cJ behavior at high frequencies either. The loss modulus calculated 
from the intra molecular stress component, as well as the direct Rouse determination 
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follow closely the square root behavior, but the contribution of "non polymer" short time 
scale inter chain forces change this behavior to almost completely mask the "polymeric" 
cross over. 
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FIG. 10: Storage modulus divided by the chain density times ksT as a function of the reduced 
frequency ujtr for the two systems. A cross over in the behavior is visible for to = 1/r/j. 
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FIG. 11: Loss modulus divided by the chain density times ksT as a function of the reduced 
frequency ojtr for the two systems (A^ = 20 - left, = 10 - right). NEMD full stress G" - 
squares, NEMD intra molecular stress G" - triangles, corrected Rouse calculation - stars, Rouse 
calculation without short times correction - small circles. A crossover in the intra molecular 
component can be seen for oj = 1/r/j. "Non polymer" relaxation mask the slope change for the 
full stress G"{uj). 



Following the discussion of the linear properties of the melt we can further explore 
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the non linear behavior studied by NEMD, considering the results presented in fig. El 
Both the relative decrease in moduli and the harmonics intensity show that the storage 
modulus has a stronger non linear behavior compared to the loss modulus at a given 
shear rate. As we have shown that G'{uj) has, at all considered frequencies, a large intra 
molecular contribution, and given that the bond potential is much steeper than the pair 
potential acting between all monomers, it is reasonable to expect that non linear effects 
due to large deformations would be more pronounced and would appear earher for the 
storage modulus. From fig. Elone can see that for u > I/tr, the higher the frequency, the 
higher the shear rate determining the onset of non-linear behavior. This critical value of 
the shear rate was found to vary linearly with frequency (701-1-') c oc uj. The system starts 
behaving as viscoelastic, and the linear regime extends to higher shear rates. In fact, as 
global chain relaxation does not take place on the time scale of the oscillations, there are 
already "frozen" slow vibration modes in the linear regime at low shear rates. Thus the 
relevant shear rate for the onset of non- linearities is shifted from l/r/j = 1/ri to 1/Tp, 
where p > 1 is the number of a higher vibration mode that still relaxes on the oscillations 
time scale at the given frequency. 

Finally, we can summarize the overall mechanical behavior of the studied polymer melt 
exhibiting several distinct regimes as shown in fig. El At low frequencies < uj < 
and low shear rate amplitudes 701^ < I/th the melt has Newtonian behavior. The viscosity 
is independent of shear rate and the elastic response is rather small as most of the stress is 
relaxed by the chain conformations. At low frequencies and high shear rates 'Jquj > l/r^f 
the system is non linear: shear thinning in viscosity, softening in moduli and harmonics in 
the measured stress. When we shift to high frequencies u > the melt exhibits more 
pronounced viscoelastic behavior with increasing storage modulus due to "frozen" modes 
and less intra molecular stress relaxation. The non linear boundary becomes frequency 
dependent and is shifted to higher shear rates determined by the relaxation time scale 
of chain segments of length < A^. The melt is then expected to exhibit glassy behavior 
at very high frequencies when no subchain relaxation whatsoever can occur within an 
oscillation. 



21 











Glassy 




Viscoelastic 


















Newtonian 


Non Linear 



FIG. 12: Mechanical behavior of the melt for different frequencies and shear rates 

V. DISCUSSION AND CONCLUSIONS 

We have discussed the visco-elastic response of a model unentangled polymer melt to 
external shear strain. We use non equilibrium molecular dynamics methods to directly 
measure the elastic moduli and the viscosity of the melt. Our aim was also to investigate 
an equilibrium based method for these quantities, inspired by Green-Kubo relations and 
possibly offering greater precision. The method we proposed is inspired by the Rouse 
model, and based on a measurement of the vibration modes of the chains at equilibrium. 
This single chain quantity can be rather accurately measured from an equilibrium simula- 
tion. Following the philosophy of the Rouse model, we use the Rouse modes of the chains 
to estimate the long times mechanical behavior of the melt and the equilibrium stress 
correlation function to estimate small times, independent of chain length non-polymeric 
behavior. The resulting model provides a satisfactory description of the viscosity and of 
the elastic moduli of the melt. The part of the different contributions in the mechanical 
stress played in the visco-elastic melt behavior was also discussed. The NEMD results 
show that inter chain interactions in the melt add a non negligible background part to the 
storage modulus and dictate the loss modulus behavior at high frequencies. These inter- 
actions should be taken into account for a precise description of the mechanical properties 
of a polymer melt made of relatively short chains. More generally, these contributions 
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are important for longer chains in several far from equilibrium situations including relax- 
ation. We measured the chain Rouse modes mean square equilibrium values and mode 
correlations directly from our simulation, so that these quantities already contain, in a 
"mean field" way, some information about the background environment of the polymer 
chains. After completing this measurements quantitatively with short times estimate of 
the stress behavior, the obtained values for the mechanical properties of the melt are 
accurate, compared to non equilibrium "direct" measurements. The NEMD results show 
that the storage modulus depends on mode relaxation over a large frequency range, the 
energy storage takes place in chain conformations whereas for the loss modulus vibra- 
tion modes are "frozen" one by one as the frequency grows higher than the inverse chain 
relaxation time and dissipation is dictated by short time scales "non polymer" stress re- 
laxation. A given vibration mode p relaxing over Tp, interpreted as the relaxation of a 
sub-chain of length N/p monomers, can contribute mostly to the loss modulus {u < l/rp) 
or to the storage modulus, if > l/rp. This picture provides an explanation for the 
overall behavior of the elastic moduli in the studied frequency range. Our results are in 
quantitative agreement with previous studies, ^ but here we focused on lower frequencies 
for the elastic moduli, based on the chain relaxation time that we estimated, in order to 
interpret the microscopic mechanisms involved. Different stress contributions on different 
time scales were measured from the simulations, thus allowing the determination of me- 
chanical response from equilibrium properties involving Rouse modes measurements and 
short times corrections. 

We examined the onset of non linear effects in the measured quantities, manifested by 
shear thinning, moduli softening and harmonics in the stress time series. The non linear 
regime is dictated by the shear rate of the solicitation and takes place at 700; > 1/tr 
for uj < 1/tr. At higher frequencies, onset of nonlinear effects is related to the strain 
amplitude rather than rate, as full chain relaxation does not take place on the time scale 
of the oscillations. Our study allows a comprehensive description of the melt mechanical 
behavior in the form of a schematic frequency - shear rate diagram shown in fig. ^| 

We finally mention that the general method presented here is not, in principle, limited 
to unentangled melts. Indeed, in the general reptation picture, the formula used for the 
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stress tensor is the same as that used in the Rouse model, formula El The relaxation of 
the modes will be dramatically slowed down by entanglement effects, so that the viscosity 
will increase rapidly. Prefactors, however, are associated with equilibrium correlation 
functions and are not affected by entanglement effects. 
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